
import pandas as pd
import numpy as np
import geopandas as gpd # Vector geospatial operations
import contextily as ctx # Used for contextual basemaps
import matplotlib.pyplot as plt
from geocube.api.core import make_geocube # Used for rasterizing
import os
import shapely
import imageio # Used for making animated GIFs
from IPython.display import Image
from osgeo import gdal # Raster operations
import zipfile
import rasterio
import rasterio.merge
import rasterio.plot
import rasterio.warp
plt.rcParams['figure.figsize'] = (20, 20)
os.listdir("input")

/usr/local/lib/python3.8/dist-packages/geopandas/_compat.py:106: UserWarning: The Shapely GEOS version (3.8.0-CAPI-1.13.1 ) is incompatible with the GEOS version PyGEOS was compiled with (3.9.0-CAPI-1.16.2). Conversions between both will be slow.
  warnings.warn(

['lds-nz-road-centrelines-topo-150k-FGDB.zip',
 'subnational-population-projections-2018base-2048.xlsx',
 'lds-nz-8m-digital-elevation-model-2012-GTiff-auckland-region.zip',
 'statsnzterritorial-authority-local-board-2021-clipped-generalised-FGDB.zip',
 'statsnzregional-council-2021-clipped-generalised-FGDB.zip',
 'lris-lcdb-v50-land-cover-database-version-50-mainland-new-zealand-FGDB.zip']


%%time
REGC = gpd.read_file("input/statsnzregional-council-2021-clipped-generalised-FGDB.zip!regional-council-2021-clipped-generalised.gdb")
AKL = REGC[REGC.REGC2021_V1_00_NAME == "Auckland Region"].copy()
# Filter out islands
AKL["geometry"] = max(AKL.geometry.explode(), key=lambda a: a.area)
# Coordinate reference system (projection)
print(AKL.crs)
# Simplify geometry to speed up clip operations
AKL = AKL.simplify(1000).buffer(1000)
ax = AKL.to_crs(epsg=3857).boundary.plot()
ax.set_title("Auckland Region clip extent")
ctx.add_basemap(ax)

epsg:2193
CPU times: user 1.84 s, sys: 878 ms, total: 2.72 s
Wall time: 14.1 s


%%time
df = gpd.read_file("zip://input/lris-lcdb-v50-land-cover-database-version-50-mainland-new-zealand-FGDB.zip!lcdb-v50-land-cover-database-version-50-mainland-new-zealand.gdb")

CPU times: user 1min 47s, sys: 2.71 s, total: 1min 50s
Wall time: 1min 50s


print(df.columns)
print(df.crs)
display(df.sample(5))

Index(['Name_2018', 'Name_2012', 'Name_2008', 'Name_2001', 'Name_1996',
       'Class_2018', 'Class_2012', 'Class_2008', 'Class_2001', 'Class_1996',
       'Wetland_18', 'Wetland_12', 'Wetland_08', 'Wetland_01', 'Wetland_96',
       'Onshore_18', 'Onshore_12', 'Onshore_08', 'Onshore_01', 'Onshore_96',
       'EditAuthor', 'EditDate', 'LCDB_UID', 'geometry'],
      dtype='object')
epsg:2193


%%time
df = gpd.clip(df, AKL)

CPU times: user 46.4 s, sys: 28.3 ms, total: 46.5 s
Wall time: 46.5 s


df.sample(5)


df.Name_2018.value_counts()

Exotic Forest                                3981
Indigenous Forest                            3673
Manuka and/or Kanuka                         2282
Broadleaved Indigenous Hardwoods             1788
Built-up Area (settlement)                   1350
High Producing Exotic Grassland              1326
Mangrove                                     1151
Urban Parkland/Open Space                    1099
Estuarine Open Water                          441
Orchard, Vineyard or Other Perennial Crop     436
Short-rotation Cropland                       362
Lake or Pond                                  326
Herbaceous Saline Vegetation                  303
Low Producing Grassland                       291
Gorse and/or Broom                            287
Forest - Harvested                            266
Sand or Gravel                                252
Deciduous Hardwoods                           201
Surface Mine or Dump                          132
Mixed Exotic Shrubland                        120
Herbaceous Freshwater Vegetation              118
Transport Infrastructure                      107
River                                          15
Flaxland                                        9
Gravel or Rock                                  9
Fernland                                        2
Matagouri or Grey Scrub                         1
Name: Name_2018, dtype: int64


def simplify_classes(code):
    if code in [1, 2, 5]:
        return 1, "Urban"
    elif code in [68,69,71]:
        return 2, "Vegetation"
    elif code in [0,20,21,22,45,46]:
        return 3, "Water"
    else:
        return 4, "Other"

summary = []
years = [1996, 2001, 2008, 2012, 2018]
for year in years:
    print(year)
    class_year = f"Class_{year}"
    df[class_year + "_simplified_code"] = df[class_year].apply(lambda c: simplify_classes(c)[0])
    df[class_year + "_simplified_name"] = df[class_year].apply(lambda c: simplify_classes(c)[1])
    summary.append(df[class_year + "_simplified_name"].value_counts())

1996
2001
2008
2012
2018


pd.DataFrame(summary).plot.area()

<AxesSubplot:>


%%capture
# %%capture suppresses output
if not os.path.isfile("land_use.gif"):
    ims = []
    years = [1996, 2001, 2008, 2012, 2018]
    for year in years:
        ax = df.plot(column=f'Class_{year}_simplified_name', legend=True)
        ax.set_title(year)
        ax.figure.tight_layout()
        canvas = ax.figure.canvas
        canvas.draw() # draw the canvas, cache the renderer
        image = np.frombuffer(canvas.tostring_rgb(), dtype='uint8')
        image = image.reshape(canvas.get_width_height()[::-1] + (3,))
        ims.append(image)
    imageio.mimsave("land_use.gif", ims, fps=1)


with open('land_use.gif','rb') as file:
    display(Image(file.read()))


cols = [f"Class_{year}_simplified_code" for year in years]
cols

['Class_1996_simplified_code',
 'Class_2001_simplified_code',
 'Class_2008_simplified_code',
 'Class_2012_simplified_code',
 'Class_2018_simplified_code']


%%time
geocube = make_geocube(
    vector_data=df,
    output_crs="epsg:2193",
    measurements=cols,
    resolution=(-100, 100),
    fill=0, # NaNs, like offshore areas, will be 0
)
geocube

CPU times: user 18.3 s, sys: 9.77 ms, total: 18.3 s
Wall time: 18.3 s

<xarray.Dataset>
Dimensions:                     (x: 1001, y: 1320)
Coordinates:
  * y                           (y) float64 6.002e+06 6.002e+06 ... 5.87e+06
  * x                           (x) float64 1.704e+06 1.704e+06 ... 1.804e+06
    spatial_ref                 int64 0
Data variables:
    Class_1996_simplified_code  (y, x) float64 0.0 0.0 0.0 0.0 ... 0.0 0.0 0.0
    Class_2001_simplified_code  (y, x) float64 0.0 0.0 0.0 0.0 ... 0.0 0.0 0.0
    Class_2008_simplified_code  (y, x) float64 0.0 0.0 0.0 0.0 ... 0.0 0.0 0.0
    Class_2012_simplified_code  (y, x) float64 0.0 0.0 0.0 0.0 ... 0.0 0.0 0.0
    Class_2018_simplified_code  (y, x) float64 0.0 0.0 0.0 0.0 ... 0.0 0.0 0.0
Attributes:
    grid_mapping:  spatial_ref

array([6002350., 6002250., 6002150., ..., 5870650., 5870550., 5870450.])

array([1703950., 1704050., 1704150., ..., 1803750., 1803850., 1803950.])

array(0)

array([[0., 0., 0., ..., 0., 0., 0.],
       [0., 0., 0., ..., 0., 0., 0.],
       [0., 0., 0., ..., 0., 0., 0.],
       ...,
       [0., 0., 0., ..., 0., 0., 0.],
       [0., 0., 0., ..., 0., 0., 0.],
       [0., 0., 0., ..., 0., 0., 0.]])

array([[0., 0., 0., ..., 0., 0., 0.],
       [0., 0., 0., ..., 0., 0., 0.],
       [0., 0., 0., ..., 0., 0., 0.],
       ...,
       [0., 0., 0., ..., 0., 0., 0.],
       [0., 0., 0., ..., 0., 0., 0.],
       [0., 0., 0., ..., 0., 0., 0.]])

array([[0., 0., 0., ..., 0., 0., 0.],
       [0., 0., 0., ..., 0., 0., 0.],
       [0., 0., 0., ..., 0., 0., 0.],
       ...,
       [0., 0., 0., ..., 0., 0., 0.],
       [0., 0., 0., ..., 0., 0., 0.],
       [0., 0., 0., ..., 0., 0., 0.]])

array([[0., 0., 0., ..., 0., 0., 0.],
       [0., 0., 0., ..., 0., 0., 0.],
       [0., 0., 0., ..., 0., 0., 0.],
       ...,
       [0., 0., 0., ..., 0., 0., 0.],
       [0., 0., 0., ..., 0., 0., 0.],
       [0., 0., 0., ..., 0., 0., 0.]])

array([[0., 0., 0., ..., 0., 0., 0.],
       [0., 0., 0., ..., 0., 0., 0.],
       [0., 0., 0., ..., 0., 0., 0.],
       ...,
       [0., 0., 0., ..., 0., 0., 0.],
       [0., 0., 0., ..., 0., 0., 0.],
       [0., 0., 0., ..., 0., 0., 0.]])


geocube.Class_2018_simplified_code.plot()

<matplotlib.collections.QuadMesh at 0x7faeb6597b20>


for year in years:
    print(year)
    outfile = f"output/land_use_{year}.tif"
    if not os.path.isfile(outfile):
        geocube[f"Class_{year}_simplified_code"].rio.to_raster(outfile, dtype=np.byte) # Use np.byte for smaller output filesize

1996
2001
2008
2012
2018


%%time
TALB = gpd.read_file("input/statsnzterritorial-authority-local-board-2021-clipped-generalised-FGDB.zip!territorial-authority-local-board-2021-clipped-generalised.gdb")
TALB = gpd.clip(TALB, AKL)
TALB.head()

CPU times: user 1.2 s, sys: 9.79 ms, total: 1.21 s
Wall time: 1.21 s


# From https://www.stats.govt.nz/information-releases/subnational-population-projections-2018base2048#map
pop = pd.concat(pd.read_excel(
    "input/subnational-population-projections-2018base-2048.xlsx",
    sheet_name=["Table 5", "Table 6"],
    skiprows=5,
    usecols="A,B,G",
    names=["area", "year", "population"],
    nrows=231
)).fillna(method='ffill').reset_index()
pop.replace("Maungakiekie-Tamaki local board area", "Maungakiekie-Tāmaki local board area", inplace=True)
pop.year = pop.year.astype(str)
pop


pop[pop.area.str.contains("local")].pivot(index='year', columns='area', values='population').plot()

<AxesSubplot:xlabel='year'>


pop = pop.pivot(index='area', columns='year', values='population')
pop


cols = pop.columns.tolist()
cols

['1996',
 '2001',
 '2006',
 '2013',
 '2018',
 '2023',
 '2028',
 '2033',
 '2038',
 '2043',
 '2048']


pop = pd.merge(TALB, pop, left_on=TALB.TALB2021_V1_00_NAME.str.lower().str.strip(), right_on=pop.index.str.lower().str.strip(), how="left")
pop[["TALB2021_V1_00_NAME"] + cols]


d = pop[cols].to_numpy()
d.min(), d.mean(), d.max()

(17800, 84495.2479338843, 190900)


%%capture
# %%capture suppresses output
if not os.path.isfile("pop.gif"):
    ims = []
    for col in cols:
        ax = pop.plot(column=col, legend=True, vmin=0, vmax=d.max(), cmap='OrRd')
        ax.set_title(col)
        ax.figure.tight_layout()
        canvas = ax.figure.canvas
        canvas.draw() # draw the canvas, cache the renderer
        image = np.frombuffer(canvas.tostring_rgb(), dtype='uint8')
        image = image.reshape(canvas.get_width_height()[::-1] + (3,))
        ims.append(image)
    imageio.mimsave("pop.gif", ims, fps=1)


with open('pop.gif','rb') as file:
    display(Image(file.read()))


%%time
popcube = make_geocube(
    vector_data=pop,
    measurements=cols,
    like=geocube, # Ensures dimensions match
    fill=0 # NaNs, like offshore areas, will be 0
)
popcube

CPU times: user 1.88 s, sys: 0 ns, total: 1.88 s
Wall time: 1.88 s

<xarray.Dataset>
Dimensions:      (x: 1001, y: 1320)
Coordinates:
  * y            (y) float64 6.002e+06 6.002e+06 ... 5.871e+06 5.87e+06
  * x            (x) float64 1.704e+06 1.704e+06 ... 1.804e+06 1.804e+06
    spatial_ref  int64 0
Data variables:
    1996         (y, x) float64 0.0 0.0 0.0 0.0 0.0 0.0 ... 0.0 0.0 0.0 0.0 0.0
    2001         (y, x) float64 0.0 0.0 0.0 0.0 0.0 0.0 ... 0.0 0.0 0.0 0.0 0.0
    2006         (y, x) float64 0.0 0.0 0.0 0.0 0.0 0.0 ... 0.0 0.0 0.0 0.0 0.0
    2013         (y, x) float64 0.0 0.0 0.0 0.0 0.0 0.0 ... 0.0 0.0 0.0 0.0 0.0
    2018         (y, x) float64 0.0 0.0 0.0 0.0 0.0 0.0 ... 0.0 0.0 0.0 0.0 0.0
    2023         (y, x) float64 0.0 0.0 0.0 0.0 0.0 0.0 ... 0.0 0.0 0.0 0.0 0.0
    2028         (y, x) float64 0.0 0.0 0.0 0.0 0.0 0.0 ... 0.0 0.0 0.0 0.0 0.0
    2033         (y, x) float64 0.0 0.0 0.0 0.0 0.0 0.0 ... 0.0 0.0 0.0 0.0 0.0
    2038         (y, x) float64 0.0 0.0 0.0 0.0 0.0 0.0 ... 0.0 0.0 0.0 0.0 0.0
    2043         (y, x) float64 0.0 0.0 0.0 0.0 0.0 0.0 ... 0.0 0.0 0.0 0.0 0.0
    2048         (y, x) float64 0.0 0.0 0.0 0.0 0.0 0.0 ... 0.0 0.0 0.0 0.0 0.0
Attributes:
    grid_mapping:  spatial_ref

array([6002350., 6002250., 6002150., ..., 5870650., 5870550., 5870450.])

array([1703950., 1704050., 1704150., ..., 1803750., 1803850., 1803950.])

array(0)

array([[0., 0., 0., ..., 0., 0., 0.],
       [0., 0., 0., ..., 0., 0., 0.],
       [0., 0., 0., ..., 0., 0., 0.],
       ...,
       [0., 0., 0., ..., 0., 0., 0.],
       [0., 0., 0., ..., 0., 0., 0.],
       [0., 0., 0., ..., 0., 0., 0.]])

array([[0., 0., 0., ..., 0., 0., 0.],
       [0., 0., 0., ..., 0., 0., 0.],
       [0., 0., 0., ..., 0., 0., 0.],
       ...,
       [0., 0., 0., ..., 0., 0., 0.],
       [0., 0., 0., ..., 0., 0., 0.],
       [0., 0., 0., ..., 0., 0., 0.]])

array([[0., 0., 0., ..., 0., 0., 0.],
       [0., 0., 0., ..., 0., 0., 0.],
       [0., 0., 0., ..., 0., 0., 0.],
       ...,
       [0., 0., 0., ..., 0., 0., 0.],
       [0., 0., 0., ..., 0., 0., 0.],
       [0., 0., 0., ..., 0., 0., 0.]])

array([[0., 0., 0., ..., 0., 0., 0.],
       [0., 0., 0., ..., 0., 0., 0.],
       [0., 0., 0., ..., 0., 0., 0.],
       ...,
       [0., 0., 0., ..., 0., 0., 0.],
       [0., 0., 0., ..., 0., 0., 0.],
       [0., 0., 0., ..., 0., 0., 0.]])

array([[0., 0., 0., ..., 0., 0., 0.],
       [0., 0., 0., ..., 0., 0., 0.],
       [0., 0., 0., ..., 0., 0., 0.],
       ...,
       [0., 0., 0., ..., 0., 0., 0.],
       [0., 0., 0., ..., 0., 0., 0.],
       [0., 0., 0., ..., 0., 0., 0.]])


for col in cols:
    outfile = f"output/pop_{col}.tif"
    if not os.path.isfile(outfile):
        # byte max value is 255, and we have larger values than that here. uint32 max value is 4294967295
        popcube[col].rio.to_raster(outfile, dtype=np.uint32)


%%time
roads = gpd.read_file("input/lds-nz-road-centrelines-topo-150k-FGDB.zip!nz-road-centrelines-topo-150k.gdb")

CPU times: user 11.8 s, sys: 10.8 ms, total: 11.8 s
Wall time: 11.8 s


%%time
akl_roads = gpd.clip(roads, AKL)

CPU times: user 23.1 s, sys: 0 ns, total: 23.1 s
Wall time: 23.1 s


# If a road has a highway number (hway_num not None), it's a highway/motorway
mway = akl_roads[~akl_roads.hway_num.isna()].copy()
mway


mway.name.value_counts().head(50).plot(kind="barh").invert_yaxis()


mway.hway_num.value_counts().head(50).plot(kind="barh").invert_yaxis()


ax = mway.to_crs(epsg=3857).plot()
ax.set_title("Auckland Region motorways")
ctx.add_basemap(ax)


%%time
mway_cube = make_geocube(
    vector_data=mway,
    measurements=["lane_count"],
    like=geocube, # Ensures dimensions match
    fill=0, # 0 works fine here, as every mway has at least one lane
)
mway_cube

CPU times: user 203 ms, sys: 162 µs, total: 203 ms
Wall time: 200 ms

<xarray.Dataset>
Dimensions:      (x: 1001, y: 1320)
Coordinates:
  * y            (y) float64 6.002e+06 6.002e+06 ... 5.871e+06 5.87e+06
  * x            (x) float64 1.704e+06 1.704e+06 ... 1.804e+06 1.804e+06
    spatial_ref  int64 0
Data variables:
    lane_count   (y, x) float64 0.0 0.0 0.0 0.0 0.0 0.0 ... 0.0 0.0 0.0 0.0 0.0
Attributes:
    grid_mapping:  spatial_ref

array([6002350., 6002250., 6002150., ..., 5870650., 5870550., 5870450.])

array([1703950., 1704050., 1704150., ..., 1803750., 1803850., 1803950.])

array(0)

array([[0., 0., 0., ..., 0., 0., 0.],
       [0., 0., 0., ..., 0., 0., 0.],
       [0., 0., 0., ..., 0., 0., 0.],
       ...,
       [0., 0., 0., ..., 0., 0., 0.],
       [0., 0., 0., ..., 0., 0., 0.],
       [0., 0., 0., ..., 0., 0., 0.]])


mway_cube.lane_count.plot()
outfile = "output/mway.tif"
if not os.path.isfile(outfile):
    mway_cube.lane_count.rio.to_raster(outfile, dtype=np.byte)


src_ds = gdal.Open("output/mway.tif")
srcband = src_ds.GetRasterBand(1)
dst_filename = "output/mway_dist.tif"
drv = gdal.GetDriverByName('GTiff')
dst_ds = drv.Create( dst_filename,
                     src_ds.RasterXSize, src_ds.RasterYSize, 1,
                     gdal.GetDataTypeByName('UInt16'))
dst_ds.SetGeoTransform( src_ds.GetGeoTransform() )
dst_ds.SetProjection( src_ds.GetProjectionRef() )
dstband = dst_ds.GetRasterBand(1)
prox = gdal.ComputeProximity(srcband,dstband,["DISTUNITS=GEO"]) # Encoded value is distance from motorway in meters
# Garbage collection of this variable flushes write
dst_ds = None
dst_ds = gdal.Open(dst_filename)
mway_dist = np.array(dst_ds.GetRasterBand(1).ReadAsArray())
print(mway_dist.shape)
plt.imshow(mway_dist)
plt.title("Distance from motorways in Auckland")
cb = plt.colorbar()
cb.ax.set_title("Distance (m)")

(1320, 1001)

Text(0.5, 1.0, 'Distance (m)')


%%time
skytower = [-36.84838748948485, 174.7621736911587]
skytower = gpd.points_from_xy(x=[skytower[1]], y=[skytower[0]])
cbd = gpd.GeoDataFrame([{"name": "Skytower", "value": 1}], geometry=skytower, crs="EPSG:4326").to_crs(epsg=2193)
display(cbd)
cbd.geometry = cbd.geometry.buffer(1000)
cbd_cube = make_geocube(
    vector_data=cbd,
    like=geocube, # Ensures dimensions match
    fill=0
)
display(cbd_cube)
outfile = "output/cbd.tif"
if not os.path.isfile(outfile):
    cbd_cube.value.rio.to_raster(outfile, dtype=np.byte)
    
src_ds = gdal.Open("output/cbd.tif")
srcband = src_ds.GetRasterBand(1)
dst_filename = "output/cbd_dist.tif"
drv = gdal.GetDriverByName('GTiff')
dst_ds = drv.Create( dst_filename,
                     src_ds.RasterXSize, src_ds.RasterYSize, 1,
                     gdal.GetDataTypeByName('UInt16'))
dst_ds.SetGeoTransform( src_ds.GetGeoTransform() )
dst_ds.SetProjection( src_ds.GetProjectionRef() )
dstband = dst_ds.GetRasterBand(1)
prox = gdal.ComputeProximity(srcband,dstband,["DISTUNITS=GEO"]) # Encoded value is distance from motorway in meters
# Garbage collection of this variable flushes write
dst_ds = None
dst_ds = gdal.Open(dst_filename)
cbd_dist = np.array(dst_ds.GetRasterBand(1).ReadAsArray())
print(cbd_dist.shape)
plt.imshow(cbd_dist)
plt.title("Distance from CBD in Auckland")
cb = plt.colorbar()
cb.ax.set_title("Distance (m)")

<xarray.Dataset>
Dimensions:      (x: 1001, y: 1320)
Coordinates:
  * y            (y) float64 6.002e+06 6.002e+06 ... 5.871e+06 5.87e+06
  * x            (x) float64 1.704e+06 1.704e+06 ... 1.804e+06 1.804e+06
    spatial_ref  int64 0
Data variables:
    value        (y, x) float64 0.0 0.0 0.0 0.0 0.0 0.0 ... 0.0 0.0 0.0 0.0 0.0
Attributes:
    grid_mapping:  spatial_ref

array([6002350., 6002250., 6002150., ..., 5870650., 5870550., 5870450.])

array([1703950., 1704050., 1704150., ..., 1803750., 1803850., 1803950.])

array(0)

array([[0., 0., 0., ..., 0., 0., 0.],
       [0., 0., 0., ..., 0., 0., 0.],
       [0., 0., 0., ..., 0., 0., 0.],
       ...,
       [0., 0., 0., ..., 0., 0., 0.],
       [0., 0., 0., ..., 0., 0., 0.],
       [0., 0., 0., ..., 0., 0., 0.]])

(1320, 1001)
CPU times: user 360 ms, sys: 10 ms, total: 370 ms
Wall time: 365 ms

Text(0.5, 1.0, 'Distance (m)')


bounds = AKL.total_bounds.tolist()
bounds

[1703081.9789640256, 5870396.320936217, 1804839.668875325, 6002367.198185163]


zf = zipfile.ZipFile('input/lds-nz-8m-digital-elevation-model-2012-GTiff-auckland-region.zip')
tiles = [file for file in zf.namelist() if file.endswith(".tif")]
tiles

['EJ.tif',
 'DM.tif',
 'EL.tif',
 'DL.tif',
 'DJ.tif',
 'FK.tif',
 'DK.tif',
 'EK.tif',
 'FL.tif']


tile_datasets = [rasterio.open(f'zip://input/lds-nz-8m-digital-elevation-model-2012-GTiff-auckland-region.zip!{tile}') for tile in tiles]
DEM, transform = rasterio.merge.merge(tile_datasets, bounds = bounds, res = (100,100), dtype=np.int16)


print(np.nanmin(DEM), np.nanmean(DEM), np.nanmax(DEM), DEM.shape)
rasterio.plot.show(np.where(DEM>=0, DEM, np.nan), cmap='terrain')

-32767 -18316.517014943143 697 (1, 1320, 1018)

<AxesSubplot:>


cbd_dist = rasterio.open("output/cbd_dist.tif")
cbd_dist.read().shape

(1, 1320, 1001)


NODATA = -32767
warped_DEM = np.full_like(cbd_dist.read(), NODATA, dtype=np.int16)

rasterio.warp.reproject(DEM,
                        warped_DEM,
                        src_transform=transform,
                        src_nodata=NODATA,
                        dst_nodata=NODATA,
                        src_crs=tile_datasets[0].crs,
                        dst_crs=cbd_dist.crs,
                        dst_transform=cbd_dist.transform
)
print(DEM.shape, warped_DEM.shape)

(1, 1320, 1018) (1, 1320, 1001)


meta = tile_datasets[0].meta
print(meta)
meta.update({
    "dtype": "int16",
    "height": cbd_dist.height,
    "width": cbd_dist.width,
    "transform": cbd_dist.transform
})
print(meta)
outfile = "output/dem.tif"
if not os.path.isfile(outfile):
    with rasterio.open(outfile, "w", **meta) as dest:
       dest.write(warped_DEM)

{'driver': 'GTiff', 'dtype': 'float32', 'nodata': -32767.0, 'width': 3028, 'height': 8192, 'count': 1, 'crs': CRS.from_epsg(2193), 'transform': Affine(8.0, 0.0, 1679712.0,
       0.0, -8.0, 5963776.0)}
{'driver': 'GTiff', 'dtype': 'int16', 'nodata': -32767.0, 'width': 1001, 'height': 1320, 'count': 1, 'crs': CRS.from_epsg(2193), 'transform': Affine(100.0, 0.0, 1703900.0,
       0.0, -100.0, 6002400.0)}


gdal.DEMProcessing('output/slope.tif', outfile, 'slope')
slope = rasterio.open("output/slope.tif").read()
print(np.nanmin(slope), np.nanmean(slope), np.nanmax(slope), slope.shape)
rasterio.plot.show(np.where(slope>=0, slope, np.nan), cmap='terrain')

-9999.0 -5778.3354 52.492313 (1, 1320, 1001)

<AxesSubplot:>


rasters = [rasterio.open(f"output/{f}") for f in os.listdir("output") if f.endswith(".tif")]
params = set([(r.shape, r.res, r.crs, r.count, r.bounds) for r in rasters])
print(params)
# Assert all rasters have the same shape, pixel size, CRS, number of bands, and bounds
assert len(params) == 1

{((1320, 1001), (100.0, 100.0), CRS.from_epsg(2193), 1, BoundingBox(left=1703900.0, bottom=5870400.0, right=1804000.0, top=6002400.0))}


!ls -Ggh output

total 66M
-rw-r--r-- 1 1.3M Apr 19 14:35 cbd.tif
-rw-r--r-- 1 2.6M Apr 19 14:35 cbd_dist.tif
-rw-r--r-- 1 2.6M Apr 19 11:42 dem.tif
-rw-r--r-- 1 1.3M Apr 15 14:25 land_use_1996.tif
-rw-r--r-- 1 1.3M Apr 15 14:25 land_use_2001.tif
-rw-r--r-- 1 1.3M Apr 15 14:25 land_use_2008.tif
-rw-r--r-- 1 1.3M Apr 15 14:25 land_use_2012.tif
-rw-r--r-- 1 1.3M Apr 15 14:25 land_use_2018.tif
-rw-r--r-- 1 1.3M Apr 15 14:37 mway.tif
-rw-r--r-- 1 2.6M Apr 19 14:18 mway_dist.tif
-rw-r--r-- 1 5.1M Apr 19 14:05 pop_1996.tif
-rw-r--r-- 1 5.1M Apr 19 14:10 pop_2001.tif
-rw-r--r-- 1 5.1M Apr 19 14:11 pop_2006.tif
-rw-r--r-- 1 5.1M Apr 19 14:11 pop_2013.tif
-rw-r--r-- 1 5.1M Apr 19 14:11 pop_2018.tif
-rw-r--r-- 1 5.1M Apr 19 14:11 pop_2023.tif
-rw-r--r-- 1 5.1M Apr 19 14:11 pop_2028.tif
-rw-r--r-- 1 5.1M Apr 19 14:11 pop_2033.tif
-rw-r--r-- 1 5.1M Apr 19 14:11 pop_2038.tif
-rw-r--r-- 1 5.1M Apr 19 14:11 pop_2043.tif
-rw-r--r-- 1 5.1M Apr 19 14:11 pop_2048.tif
-rw-r--r-- 1 5.1M Apr 19 14:18 slope.tif

	TALB2021_V1_00	TALB2021_V1_00_NAME	TALB2021_V1_00_NAME_ASCII	LAND_AREA_SQ_KM	AREA_SQ_KM	Shape_Length	geometry
2	00300	Kaipara District	Kaipara District	3108.960347	3108.960347	934117.838186	MULTIPOLYGON (((1718576.593 5980699.195, 17185...
4	01200	Hauraki District	Hauraki District	1270.133527	1270.133527	311282.487609	POLYGON ((1801726.041 5899678.947, 1803788.474...
5	01300	Waikato District	Waikato District	4403.797301	4451.122298	680627.805117	POLYGON ((1747419.052 5870470.171, 1747413.783...
44	07601	Rodney Local Board Area	Rodney Local Board Area	2275.114945	2275.114945	920471.582397	MULTIPOLYGON (((1717309.547 5972989.294, 17173...
45	07602	Hibiscus and Bays Local Board Area	Hibiscus and Bays Local Board Area	110.097301	110.097301	162589.255306	MULTIPOLYGON (((1751233.745 5948474.046, 17512...

	level_0	level_1	area	year	population
0	Table 5	0	Far North district	1996	54500
1	Table 5	1	Far North district	2001	56400
2	Table 5	2	Far North district	2006	57500
3	Table 5	3	Far North district	2013	60600
4	Table 5	4	Far North district	2018	67900
...	...	...	...	...	...
457	Table 6	226	Franklin local board area	2028	96900
458	Table 6	227	Franklin local board area	2033	109400
459	Table 6	228	Franklin local board area	2038	122000
460	Table 6	229	Franklin local board area	2043	134500
461	Table 6	230	Franklin local board area	2048	146900

year	1996	2001	2006	2013	2018	2023	2028	2033	2038	2043	2048
area
Albert-Eden local board area	86100	90000	96200	100000	103700	106900	111200	116500	121600	126300	130600
Auckland	1115800	1218300	1373000	1493200	1654800	1778700	1891800	2001800	2107000	2207800	2302900
Devonport-Takapuna local board area	51600	52300	55700	58500	60500	62700	65700	68800	71800	74800	77600
Far North district	54500	56400	57500	60600	67900	72400	75000	77300	79100	80600	81700
Franklin local board area	48400	52900	60600	68300	77700	85700	96900	109400	122000	134500	146900
Gisborne district	47200	45500	45900	47000	49500	51900	53100	54000	54600	55000	55200
Great Barrier local board area	1230	1100	940	950	960	1030	1080	1120	1150	1180	1200
Hamilton city	113500	121200	134800	150200	168600	183000	194400	205400	216000	226500	236600
Hauraki district	18550	18000	18300	18600	20700	21800	22200	22400	22300	22100	21800
Henderson-Massey local board area	82200	90900	103700	113500	124600	134100	142500	150500	158200	165700	172900
Hibiscus and Bays local board area	68000	75100	85200	94000	108500	117700	123500	126500	128900	130900	132400
Howick local board area	80900	97300	119100	135000	149400	160500	167800	174400	180400	185900	190900
Kaipara district	17800	17950	18550	20500	23700	25900	27200	28300	29100	29800	30300
Kaipātiki local board area\n	74700	78200	83600	87000	92900	97300	99400	101100	102500	103600	104700
Kawerau district	8120	7290	7150	6650	7460	7910	8000	8020	7970	7860	7720
Manurewa local board area	60100	69100	81600	87000	100900	109600	114000	118000	121700	124900	127600
Matamata-Piako district	30300	30300	31200	32900	35300	37000	37900	38700	39200	39500	39600
Maungakiekie-Tāmaki local board area	61700	66400	70400	73700	80500	86600	93800	101300	108800	116100	123000
Māngere-Ōtāhuhu local board area\n	61400	64600	72500	75300	82700	88700	92500	96300	99800	103100	105900
Papakura local board area	38000	38900	43100	48200	61100	70000	77000	81700	86100	90500	95000
Puketāpapa local board area\n	46100	49800	53900	56300	60900	64600	68800	73000	77100	81000	84600
Rodney local board area	39000	44100	51000	57300	69100	77600	87600	99600	111700	123800	135800
Rotorua district	66600	66900	68100	68400	74800	78900	80700	82200	83400	84200	84800
South Waikato district	25800	24200	23200	23200	24800	25800	26400	26800	27000	27100	27100
Taupō district\n	31600	32500	33400	34800	38600	40900	42000	42800	43300	43600	43800
Tauranga city	79900	93600	107000	119900	142100	156900	166300	175000	183300	191400	199100
Thames-Coromandel district	25400	25800	26700	27300	30700	32400	33100	33500	33500	33300	32800
Upper Harbour local board area	24300	33700	45100	56800	66800	75900	86000	95900	105800	115400	124900
Waiheke local board area	6680	7560	8150	8630	9360	9830	10250	10650	10900	11100	11200
Waikato district	52000	53700	59500	66500	78200	86100	92800	99300	105700	111900	117700
Waipa district	38400	40000	43700	48700	55000	59300	62100	64600	66900	68900	70700
Waitematā local board area\n	44900	52500	66600	81300	88500	94700	102800	111900	120700	129200	137400
Waitomo district	10000	9780	9680	9340	9630	9740	9730	9670	9540	9330	9070
Waitākere Ranges local board area\n	42100	44200	47500	50700	54200	56900	58300	59500	60300	60900	61100
Western Bay of Plenty district	35600	38900	42900	45400	53300	58100	60900	63300	65200	66700	68000
Whakatāne district\n	34200	34100	34500	34200	37100	38800	39300	39500	39500	39300	38900
Whangārei district\n	68400	70000	76500	83700	94100	101000	105500	109500	113100	116400	119300
Whau local board area	59700	66200	73100	76700	84100	89500	95000	100500	105900	111000	115900
Ōpōtiki district\n	9620	9490	9200	8780	9670	10250	10350	10400	10300	10150	9910
Ōrākei local board area\n	70900	73000	78400	83700	87700	91300	96200	101100	105800	110300	114500
Ōtara-Papatoetoe local board area\n	67800	70400	76600	80300	90500	97700	101500	104000	106000	107500	108600
Ōtorohanga district\n	9960	9590	9310	9590	10500	11050	11300	11550	11750	11900	12000

	TALB2021_V1_00_NAME	1996	2001	2006	2013	2018	2023	2028	2033	2038	2043	2048
0	Kaipara District	17800	17950	18550	20500	23700	25900	27200	28300	29100	29800	30300
1	Hauraki District	18550	18000	18300	18600	20700	21800	22200	22400	22300	22100	21800
2	Waikato District	52000	53700	59500	66500	78200	86100	92800	99300	105700	111900	117700
3	Rodney Local Board Area	39000	44100	51000	57300	69100	77600	87600	99600	111700	123800	135800
4	Hibiscus and Bays Local Board Area	68000	75100	85200	94000	108500	117700	123500	126500	128900	130900	132400
5	Upper Harbour Local Board Area	24300	33700	45100	56800	66800	75900	86000	95900	105800	115400	124900
6	Kaipātiki Local Board Area	74700	78200	83600	87000	92900	97300	99400	101100	102500	103600	104700
7	Devonport-Takapuna Local Board Area	51600	52300	55700	58500	60500	62700	65700	68800	71800	74800	77600
8	Henderson-Massey Local Board Area	82200	90900	103700	113500	124600	134100	142500	150500	158200	165700	172900
9	Waitākere Ranges Local Board Area	42100	44200	47500	50700	54200	56900	58300	59500	60300	60900	61100
10	Waitematā Local Board Area	44900	52500	66600	81300	88500	94700	102800	111900	120700	129200	137400
11	Whau Local Board Area	59700	66200	73100	76700	84100	89500	95000	100500	105900	111000	115900
12	Albert-Eden Local Board Area	86100	90000	96200	100000	103700	106900	111200	116500	121600	126300	130600
13	Puketāpapa Local Board Area	46100	49800	53900	56300	60900	64600	68800	73000	77100	81000	84600
14	Ōrākei Local Board Area	70900	73000	78400	83700	87700	91300	96200	101100	105800	110300	114500
15	Maungakiekie-Tāmaki Local Board Area	61700	66400	70400	73700	80500	86600	93800	101300	108800	116100	123000
16	Howick Local Board Area	80900	97300	119100	135000	149400	160500	167800	174400	180400	185900	190900
17	Māngere-Ōtāhuhu Local Board Area	61400	64600	72500	75300	82700	88700	92500	96300	99800	103100	105900
18	Ōtara-Papatoetoe Local Board Area	67800	70400	76600	80300	90500	97700	101500	104000	106000	107500	108600
19	Manurewa Local Board Area	60100	69100	81600	87000	100900	109600	114000	118000	121700	124900	127600
20	Papakura Local Board Area	38000	38900	43100	48200	61100	70000	77000	81700	86100	90500	95000
21	Franklin Local Board Area	48400	52900	60600	68300	77700	85700	96900	109400	122000	134500	146900

	t50_fid	name_ascii	macronated	name	hway_num	rna_sufi	lane_count	way_count	status	surface	geometry
512	100120610	KAIPARA COAST HIGHWAY	N	KAIPARA COAST HIGHWAY	16	3007739	2	None	None	sealed	LINESTRING (1732000.000 5944172.070, 1732048.5...
2933	3198057	STATE HIGHWAY 1	N	STATE HIGHWAY 1	1	3027695	2	None	None	sealed	LINESTRING (1748581.508 5968975.145, 1748558.4...
2934	3198059	STATE HIGHWAY 1	N	STATE HIGHWAY 1	1	3027695	2	None	None	sealed	LINESTRING (1748171.047 5971284.152, 1748129.9...
3320	3200754	PAERATA ROAD	N	PAERATA ROAD	22	3000260	2	None	None	sealed	LINESTRING (1767236.112 5888088.508, 1767244.3...
3324	3200792	UPPER HARBOUR MOTORWAY	N	UPPER HARBOUR MOTORWAY	18	3047073	4	None	None	sealed	LINESTRING (1747954.314 5927269.837, 1747970.0...
...	...	...	...	...	...	...	...	...	...	...	...
138240	100048291	AUCKLAND-WAIWERA MOTORWAY	N	AUCKLAND-WAIWERA MOTORWAY	1	3067966	7	None	None	sealed	LINESTRING (1755881.018 5922863.734, 1755886.4...
138301	100048432	AUCKLAND-HAMILTON MOTORWAY	N	AUCKLAND-HAMILTON MOTORWAY	1	3017109	1	None	None	sealed	LINESTRING (1765115.647 5909916.697, 1765092.7...
138337	100048532	STATE HIGHWAY 1	N	STATE HIGHWAY 1	1	3027695	4	None	None	sealed	LINESTRING (1748892.089 5949596.727, 1748892.0...
138369	100048589	PORT ALBERT ROAD	N	PORT ALBERT ROAD	16	3013274	2	None	None	sealed	LINESTRING (1734173.019 5980575.187, 1734175.6...
138680	100118365	SOUTH-WESTERN MOTORWAY	N	SOUTH-WESTERN MOTORWAY	20	3018532	4	None	None	sealed	LINESTRING (1760066.252 5908184.133, 1760043.9...

Land use¶

Population density¶

Distance from major roads¶

Distance from CBD¶

Slope¶

	Name_2018	Name_2012	Name_2008	Name_2001	Name_1996	Class_2018	Class_2012	Class_2008	Class_2001	Class_1996	...	Wetland_96	Onshore_18	Onshore_12	Onshore_08	Onshore_01	Onshore_96	EditAuthor	EditDate	LCDB_UID	geometry
313790	Landslide	Landslide	Landslide	Landslide	Landslide	12	12	12	12	12	...	no	yes	yes	yes	yes	yes	Terralink	2004-06-30T00:00:00	lcdb2000005253	MULTIPOLYGON (((1160923.451 4898204.158, 11609...
214669	Gravel or Rock	Gravel or Rock	Gravel or Rock	Gravel or Rock	Gravel or Rock	16	16	16	16	16	...	no	yes	yes	yes	yes	yes	Terralink	2004-06-30T00:00:00	lcdb1000014095	MULTIPOLYGON (((1887044.793 5620048.149, 18870...
388703	Indigenous Forest	Indigenous Forest	Indigenous Forest	Indigenous Forest	Indigenous Forest	69	69	69	69	69	...	no	yes	yes	yes	yes	yes	Landcare Research	2014-06-30T00:00:00	lcdb1000182838	MULTIPOLYGON (((1960574.060 5784716.815, 19605...
107539	Exotic Forest	Exotic Forest	Exotic Forest	Exotic Forest	Exotic Forest	71	71	71	71	71	...	no	yes	yes	yes	yes	yes	Terralink	2004-06-30T00:00:00	lcdb1000202128	MULTIPOLYGON (((1945776.018 5602220.374, 19457...
75253	Broadleaved Indigenous Hardwoods	Broadleaved Indigenous Hardwoods	Broadleaved Indigenous Hardwoods	Broadleaved Indigenous Hardwoods	Broadleaved Indigenous Hardwoods	54	54	54	54	54	...	no	yes	yes	yes	yes	yes	Landcare Research	2004-06-30T00:00:00	lcdb2000135128	MULTIPOLYGON (((1439088.933 5045590.762, 14391...

	Name_2018	Name_2012	Name_2008	Name_2001	Name_1996	Class_2018	Class_2012	Class_2008	Class_2001	Class_1996	...	Wetland_96	Onshore_18	Onshore_12	Onshore_08	Onshore_01	Onshore_96	EditAuthor	EditDate	LCDB_UID	geometry
257236	High Producing Exotic Grassland	High Producing Exotic Grassland	Exotic Forest	Exotic Forest	Exotic Forest	40	40	71	71	71	...	no	yes	yes	yes	yes	yes	Landcare Research	2014-06-30T00:00:00	lcdb1000408113	POLYGON ((1776033.445 5897655.722, 1776027.801...
429399	Exotic Forest	Exotic Forest	Exotic Forest	Exotic Forest	Exotic Forest	71	71	71	71	71	...	no	yes	yes	yes	yes	yes	Terralink	2004-06-30T00:00:00	lcdb1000217879	POLYGON ((1779737.510 5905755.113, 1779733.434...
281720	Built-up Area (settlement)	Built-up Area (settlement)	Built-up Area (settlement)	Built-up Area (settlement)	Urban Parkland/Open Space	1	1	1	1	2	...	no	yes	yes	yes	yes	yes	Landcare Research	2011-06-30T00:00:00	lcdb1000003636	POLYGON ((1746719.161 5913844.184, 1746732.995...
365670	Manuka and/or Kanuka	Manuka and/or Kanuka	Manuka and/or Kanuka	Manuka and/or Kanuka	Manuka and/or Kanuka	52	52	52	52	52	...	no	yes	yes	yes	yes	yes	Landcare Research	2015-06-30T00:00:00	lcdb1000456637	POLYGON ((1731793.016 5938194.674, 1731764.232...
296981	Exotic Forest	Exotic Forest	Exotic Forest	Exotic Forest	Exotic Forest	71	71	71	71	71	...	no	yes	yes	yes	yes	yes	Landcare Research	2014-06-30T00:00:00	lcdb1000409942	POLYGON ((1756300.691 5879431.457, 1756287.981...